arXiv: 1509.0473Ivl [cond-mat.mtrl-sci] 15 Sep 2015 


Beyond packing of hard spheres: The effects of core softness, non-additivity, 
intermediate-range repnlsion, and many-body interactions on the glass-forming ability 

of bnlk metallic glasses 

Kai Zhang*,Meng Fan*,^’^ Yanhui Jan Schroers,^’^ Mark D. Shattuck,^’^ and Corey S. O’Hern^’^’^’® 

^Department of Mechanical Engineering and Materials Science, 

Yale University, New Haven, Connecticut, 06520, USA 
^Center for Research on Interface Structures and Phenomena, 

Yale University, New Haven, Connecticut, 06520, USA 
^Department of Physics and Benjamin Levich Institute, 

The City College of the City University of New York, New York, New York, 10031, USA 
^Department of Physics, Yale University, New Haven, Connecticut, 06520, USA 
^Department of Applied Physics, Yale University, New Haven, Connecticut, 06520, USA 

(Dated: September 17, 2015) 

When a liquid is cooled well below its melting temperature at a rate that exceeds the critical 
cooling rate Rc, the crystalline state is bypassed and a metastable, amorphous glassy state forms 
instead. Rc (or the corresponding critical casting thickness dc) characterizes the glass-forming ability 
(GFA) of each material. While silica is an excellent glass-former with small Rc < 10“^ K/s, pure 
metals and most alloys are typically poor glass-formers with large Rc > 10^° K/s. Only in the past 
thirty years have bulk metallic glasses (BMGs) been identified with Rc approaching that for silica. 

Recent simulations have shown that simple, hard-sphere models are able to identify the atomic size 
ratio and number fraction regime where BMGs exist with critical cooling rates more than 13 orders 
of magnitude smaller than those for pure metals. However, there are a number of other features 
of interatomic potentials beyond hard-core interactions. How do these other features affect the 
glass-forming ability of BMGs? In this manuscript, we perform molecular dynamics simulations 
to determine how variations in the softness and non-additivity of the repulsive core and form of 
the interatomic pair potential at intermediate distances affect the GFA of binary alloys. These 
variations in the interatomic pair potential allow us to introduce geometric frustration and change 
the crystal phases that compete with glass formation. We also investigate the effect of tuning the 
strength of the many-body interactions from zero to the full embedded atom model on the GFA 
for pure metals. We then employ the full embedded atom model for binary BMGs and show that 
hard-core interactions play the dominant role in setting the GFA of alloys, while other features of 
the interatomic potential only change the GFA by one to two orders of magnitude. Despite their 
perturbative effect, understanding the detailed form of the intermetallic potential is important for 
designing BMGs with cm or greater casting thickness. 

PACS numbers: 64.70.pe,64.70.Q-,61.43.Fs,61.66.Dk,61.43.Dq 


1. INTRODUCTION 

When metallic liquids are cooled at rates R exceeding 
the critical cooling rate Rc, crystallization can be by¬ 
passed and amorphous alloys are formed [l|. Pure met¬ 
als and most alloys are extremely poor glass formers with 
Rc > 10^° K/s. In contrast, a number of bulk metallic 
glasses (BMGs) have been identified with Rc < 1 K/s and 
critical casting thickness dc > I cm, which enables them 
to be employed in commercial applications ii- The 
discovery of novel BMGs with optimized casting thick¬ 
ness and mechanical properties has largely been a trial- 
and-error process 0, . Although combinatorial deposi¬ 

tion and characterization techniques @0 now allow effi¬ 
cient exploration of parameter space, there are an expo¬ 
nentially large number of possible BMG-forming atomic 
compositions Q. Thus, a quantitative and predictive 
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understanding of the GFA of BMG-forming alloys is nec¬ 
essary to narrow down the vast parameter space. 

Silica and polymers possess critical cooling rates that 
are more than 15 and 10 orders of magnitude lower, re¬ 
spectively, than those for pure metals (Fig. [T]). Network 
bonding in silica and chain entanglement in polymers 
provide the physical mechanisms to inhibit crystalliza¬ 
tion EMI. In contrast, the main source of geomet¬ 
ric frustration in alloys is the mismatch between atomic 
sizes (Tol - f^ . Molecular dynamics simulations of binary 
hard spheres have shown that tuning the atomic size ra¬ 
tio can decrease Rc by more than 13 orders of magni¬ 
tude [l^. Packing of hard spheres can also rationalize 
the correlation between the number of components, their 
atomic size ratios, and the GFA of BMGs [^. 

Although the packing of hard spheres plays an impor¬ 
tant role in determining the GFA of alloys, it is obvi¬ 
ous that metals possess additional features that are not 
represented by hard-sphere interactions. Other features 
of metallic interactions, such as metallic bonding [^ . 
the form of the interatomic pair potential, and many- 
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FIG. 1: Schematic diagram of crystalline order (such as bond 
orientational order i) versus the cooling rate R in K/s for 
several materials. The critical cooling rate Rc at which there 
is a rapid rise in the crystalline order is inversely correlated 
with the material’s critical casting thickness d^. Smaller Rc 
(and larger dc) indicate enhanced glass-forming ability (GF.^. 
Pure metals, e.g. Ta, are extremely poor glass formers [3. 
The GFA of the first fabricated metallic alloy Au 8 oSi 20 )10l| 
is similar to that of water 0, but is a poor glass-former 
compared to polymers [0 and silica [l^ . The best bulk 
metallic glasses (BMGs), e.g. Pd4oGu3oNiioP2o [3, possess 
cm or greater critical casting thicknesses and < 1 K/s criti¬ 
cal cooling rates (solid gray bars). In recent simulations, we 
have shown that hard-core atomic interactions can account 
for more than 13 orders of magnitude variation in Rc (thick 
dashed line) from 10° K/s for typical BMGs to 10^° K/s for 
pure metals 0. 


body interactions [3: can change the crystalline struc¬ 
ture that competes with glass formation and change the 
prediction of Rc by several orders of magnitude from the 
hard-sphere value. Compared to the ~ 13 orders of mag¬ 
nitude variation in Rc that results from the packing of 
hard-spheres, changes to Rc are small, but not negligible 
and may explain the crucial differences between an amor¬ 
phous film and a bulk metallic glass. Since the critical 
casting thickness dc is negatively correlated with Rc and 
increasing Rc by two orders of magnitude can reduce dc 
by one order of magnitude [3: more accurate models 
of intermetallic potentials are needed to identify BMGs 
with dc > 1 cm (Fig. [T|). 

The interatomic potential in the embedded atom 
model (EAM) is frequently implemented in computa¬ 
tional studies of the structural and mechanical proper¬ 
ties, as well as the dynamics, of metallic systems [0. 
The EAM potential energy includes a pairwise-additive 
term, which is in general different from the hard-sphere 
and Lennard-Jones pair potentials (Fig. [2] (a)), and a 
many-body contribution from the electron charge den¬ 


sity, which is fitted to ab initio calculations of lattice 
parameters, elastic constants, and other thermodynamic 
properties (m. [^. 

In this manuscript, we seek to identify the key fea¬ 
tures of the pairwise and many-body interactions that 
strongly influence the GFA of alloys. For example, we 
investigate the effects of the softness of the pairwise re¬ 
pulsive core, pairwise non-additivity, and the form of the 
pairwise intermediate-range repulsion on the GFA. We 
then measure the GFA for the full embedded atom mod¬ 
els of several pure metals and BMGs to determine the 
contribution of the many-body interactions on the GFA. 
We find that the changes in the GFA arising from vari¬ 
ations in the pair and many-body contributions of the 
embedded atom model are small compared to the 13 or¬ 
ders of magnitude change in GFA between monoatomic 
and binary and ternary hard-sphere systems. However, 
these peturbations to the GFA may still be important for 
discovering new bulk metallic glass formers. 

The manuscript includes three additional sections af¬ 
ter the introduction. First, in Sec. [21 we describe the 
hard-sphere, repulsive Lennard-Jones, Lennard-Jones, 
and Dzugutov-Shi potentials used to vary the form and 
non-additivity of the pairwise interactions. We also in¬ 
troduce the embedded atom model for pure metals and 
alloys. For each interatomic potential, we discuss the 
methods employed to measure the critical cooling rate 
Rc. We then report the results for the GFA for all inter¬ 
action potentials in Sec. [3] We conclude the manuscript 
in Sec. 01 


2. MODELS AND METHODS 

As described above, the embedded atom model for 
metallic systems includes pairwise and many-body inter¬ 
actions. In this section, we define three metrics (core 
softness, non-additivity, and intermediate-range repul¬ 
sion) to characterize the form of the pairwise interac¬ 
tions. We describe molecular dynamics simulations of 
monodipserse and binary systems interacting via gener¬ 
alized Lennard-Jones or Dzugutov-Shi |0 potentials 
to quantify the effects of the softness of the repulsive 
core and strength of the intermediate-range repulsion on 
the GFA. We also introduce molecular dynamics simu¬ 
lations of binary hard spheres to study variations in the 
GFA from non-additive pairwise interactions. We esti¬ 
mate values for the pairwise core softness, non-additivity, 
and form of the intermediate-range repulsive interactions 
from fits to the pairwise contributions of the EAM for 
pure metals and binary BMGs. We also introduce the 
Lennard-Jones and full EAM potentials that we employ 
to study the effects of many-body interactions on the 
GFA. 
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FIG. 2: (a) The pairwise potentials u(rij) (in eV) as a function of interatomic separation Vij for Zr-Zr (solid line), Cu-Cu 
(dotted line), and Zr-Cu (dashed line) interactions for the embedded atom model for Zr-Cu alloys [2^. (b) Generalized 

Lennard-Jones (Eq. [T|) (dashed lines) and repulsive Lennard-Jones (Eq. [2j (dotted lines) interatomic potentials for several 
values of the core softness exponent m — 1, 3, 5, 8, and 12 (from left to right) compared to the hard-sphere potential (thick 
solid line), (c) Dzugutov-Shi interatomic potential (Eq. O (solid line) decomposed into the Lennard-Jones (dotted line) and 
sinusoidal “bump” potentials (dashed line). 


2.1. Lennard-Jones (LJ) and Repulsive tential (Fig. [2] (b)), 

Lennard-Jones (RLJ) Potentials 

To tune the softness of the pairwise repulsive core , 

we employ the generalized m-n Lennard-Jones (LJ) po- 
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LJ 


N) = 
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where atj = {ai + aj)/2, Ui is the diameter of atom i, and 
e is the energy scale of the interaction. The interaction 
potential has a minimum Um = at 
The exponent m (or equivalently the curvature k of the 
pair potential at the minimum) controls the softness of 
the repulsive core, where smaller m corresponds to softer 
interactions. Note that the generalized Lennard-Jones 


potential is fixed at ULj(ry ) = 


— „, 12-6 


Tj 


(ry ) for rij > 


To separate the effects of the attractive interactions from 
the repulsive core, we also studied the gene ralized m-n 
repulsive Lennard-Jones (RLJ) potential 33| as shown in 
Fig. [21(b): 


u 


m—n 

RLJ 


i^ij) 


0 ^ '^ij ^ '^m- 


( 2 ) 


To obtain physical values for the softness exponent m, 
we fit the repulsive part of the EAM pair potential of 
typical BMG-forming elements to (r). As shown in 
Table U we find that m varies from approximately 3 to 
14. The repulsive cores for most metals appear softer 
than Lennard-Jones interactions with m = 12. 

To investigate the effects of softness of the pairwise 
repulsive core on the GFA of metallic systems, we per¬ 
formed molecular dynamics (MD) simulations of = 
1372 spherical atoms with mass mo that interact via the 
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rij < Tm 
rij ^ 


( 1 ) 


generalized Lennard-Jones and repulsive Lennard-Jones 
potentials with n = 6 and a range of m values. We 
studied three binary LJ systems with softness exponents 
mA = ms = mAB = 12 (LJ12-6), m^ = ms = mAB = 5 
(LJ5-6), and m^ = 12, m^ = 5, and mAB = 8 
(LJ12-6/LJ5-6). We set the atomic diameter ratio to 
be a = (JbIcta = 0.95 and varied the number fraction 
of small atoms xb = Nb/N from 0 to 1. Temperatures 
and times are given in units of e/fcs and aA^/moje, re¬ 
spectively. After equilibrating the systems at high tem¬ 
perature To = 2, the liquids were cooled exponentially 
T(t) = Toexp(—i?t) with rate R to low temperature, 
Tf = 0.01, using the Gaussian constraint thermostat 
with time step At = 0.001. Constant volume V sim¬ 
ulations at number density pa\ = Na\/V = 1 were 
performed for both the LJ and RLJ models. For LJ sys¬ 
tems, we also cooled systems with the constraint that 
the pressure p (in units of e/cr^) decreased exponentially 
in time from an initial pressure po = 1 to final pressure 
Pf = 0.001 according to 


p{t) = po exp 


log(po/p/) ' 
log{To/Tf) _ 


( 3 ) 


using a Gaussian constraint barostat . A cooling rate 
of i? = 1 in the units used in the MD simulations corre- 
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TABLE I: Softness exponent m from the re¬ 
pulsive Lennard-Jones potential (Eq. [S} for the 
self-part of the pair potential contribution to 
the embedded atom model for common atomic 
species found in BMGs [^ . The exponent m 
varies linearly with the curvature k (given in 
units of c/o-a) of the interatomic potential at its 
minimum Vm- References for the EAM poten¬ 
tials are provided in columns 4 and 8. Several 
atom types possess multiple EAM potentials. 


atom 

K, 

m 

Ref. 

atom 

K, 

m 

Ref. 

Zr 

38.77 

8.14 

L29J 

Pb 

35.30 

7.41 

m 


31.77 

6.67 

[341 

Mg 

42.69 

8.97 

[341 


66.91 

14.05 

[281 

Fe 

31.41 

6.60 

[341 

Ag 

48.77 

10.24 

[291 

Co 

32.88 

6.90 

[341 


35.49 

7.45 

[341 


38.06 

7.99 

m 


40.35 

8.47 

[361 

Ta 

21.36 

4.49 

m 

Al 

33.48 

7.03 

[291 


15.64 

3.28 

[37,1 


20.65 

4.34 

[341 


24.56 

5.16 

[381 


10.89 

2.29 

[39| 

Cu 

28.66 

6.02 

[341 

Ni 

34.36 

7.21 

[29] 

An 

38.40 

8.06 

m 


30.66 

6.44 

[341 


48.20 

10.12 

[401 


47.96 

10.07 

m 

Ti 

32.93 

6.92 

m 

Pd 

43.72 

9.18 

[29] 

Mo 

20.42 

4.29 

m 


33.39 

7.01 

[341 

W 

22.63 

4.75 

[341 

Pt 

43.47 

9.13 

[29] 

Nb 

28.19 

5.92 

[42] 


24.04 

5.05 

m 





spends to a cooling rate of 10^® K/s using cta ~ 3 x 10“^° 
m, e/ks ~ 10^ K, and molar mass M ~ 10“^ kg/mol, 
which are typical values for BMGs [3^ . 


2.2. Non-additive binary hard spheres 

The sizes of metallic atoms are often estimated from 
the first peak of the radial distribution function g{r) of 
crystalline and disordered solids Si- In binary alloys 
with species A and B, the repulsive core aAB between 
atoms A and B can differ from the average diameter 
ctab = {o'A + Ob) 12. We quantify the non-additivity 
of the pairwise repulsive core using the parameter 

S = ^ - 1. (4) 

O AB 

Many binary alloys possess E < 0, which indicates that 
the repulsive core oab between A and B atoms is smaller 
than the average diameter. We list a a, ob, oab, and S 
for several binary alloys obtained from EAM calculations 
of g{r) in Table[TTl Non-additive binary hard spheres have 
been shown to form exotic crystalline structures, in par¬ 
ticular intermetallic compounds [i^, In addition, 

non-additivity due to bond shortening with E < 0 can 
lead to unusual intermediate-range order in BMGs [13- 
Id^ . The well-studied Kob-Andersen model for Ni 8 oP 20 
glasses also has E = —0.149 [50l |. 

To study the effects of nonadditivity on the GFA, we 
compressed N = 500 binary hard spheres with mass toq 


TABLE II: Atomic diameters (cta 
and (tb in A) determined by the first 
peak of the radial distribution func¬ 
tion g{r) obtained from EAM sim¬ 
ulations of several binary alloys [2^. 
We also list ctab from g{r), the diam¬ 
eter ratio a, and the non-additivity 
parameter E. 


Alloy 

OA 

Ob 

Oab 

a 

E 

Zr-Cu 

3.15 

2.49 

2.75 

0.79 

-0.025 

Ni-P 

2.57 

2.19 

2.23 

0.85 

-0.063 

Zr-Ni 

3.23 

2.43 

2.69 

0.75 

-0.049 

Zr-Al 

3.21 

2.69 

2.93 

0.84 

-0.007 

Ag-Al 

2.87 

2.69 

2.69 

0.94 

-0.032 

Mg-Cu 

3.11 

2.47 

2.69 

0.79 

-0.05 

Mg-Ti 

2.97 

2.77 

2.99 

0.93 

0.042 

Y-Mg 

3.51 

3.03 

3.27 

0.86 

0 

Pd-Si 

2.97 

2.39 

2.51 

0.80 

-0.063 

Zr-Pt 

3.39 

2.91 

2.73 

0.86 

-0.1333 

Cu-Ni 

2.51 

2.45 

2.49 

0.98 

0.004 

Mg-Al 

2.99 

2.81 

2.99 

0.94 

0.031 


that interact pairwise via 


UHsifij) 


: ‘^ij^ Oij 
0 ^ T^j ^ Oij 


(5) 


over a range of diameter ratios a and number fractions 
of the small sphere xb using event-driven MD simula¬ 
tions. We first equilibrated liquid states at packing frac¬ 
tion cj) = 0.25. To compress the system, we ran the 
MD simulations at constant volume for a time interval 
T, and then compressed the system instantaneously until 
the closest pair of spheres came into contact [H, . We 

performed successive compressions until the pressure in¬ 
creased to 10^, which corresponds to {4'j — (l>)/4>j < 10“^, 
where 4>j is the packing fraction at the onset of jamming. 
We varied the compression rate i? = l/r over 5 orders of 
magnitude [T^. We report R in units of ^/UbT/ mQo\. 
Note that in these units i? = 1 corresponds to a cooling 
rate of 10^^ K/s for alloys (Fi| . 


2.3. Dzugutov-Shi (DZ) potential 

The pair potential of many metallic systems includes 
intermediate-range repulsive interactions in addi¬ 
tion to short-range attractive interactions, which can 
give rise to intermediate-range positional order [sl, . 
Intermediate-range pairwise repulsive interactions are of¬ 
ten modeled using the Dzugutov potential [13, 

Shi et. al. introduced a modified version of the origi¬ 
nal Dzugutov potential that allows one to continuously 
tune the interaction potential between the LJ potential 
to one that includes intermediate-range repulsion M- 
The Dzugutov-Shi (DZ) potential is given by 

{'^ij) 4“ '^bump ) , 


( 6 ) 
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TABLE III: Values of the 
parameters A, and 5 
(Eq. that describe the 
strength and range of the 
Dzugutov-Shi interatomic 
potential fit to the self-part 
of the pair potential of 
the embedded atom model 
for several atomic species. 
The fifth column provides 
references for the EAM for 
each atom type. 


atom 

C (eV) 

A 

5 

Ref. 

Zr 

0.42 

1.16 

2.24 

[29J 

Ag 

0.16 

1.29 

2.20 

[291 

Cu 

0.43 

1.18 

1.73 

[29] 

Ni 

0.38 

1.19 

1.72 

[29] 

A1 

0.10 

1.76 

2.35 

[29] 


0.26 

1.29 

1.99 

[39] 

Co 

0.12 

1.56 

2.66 

[3^ 


where the “bump” potential Ubump(?'y) models the 
intermediate-range repulsive interactions using a sinu¬ 
soidal pulse, 


'^bump ij ) 



A < TijlOij < (5 
otherwise, 


(7) 

of the strength ^ within the range Xuij < rij < Saij. 
The location of the peak and width of Ubump are given 
by (A-|-^)/2 and S— X. To obtain physical values for A, 
and i5, we fit the DZ potential to the EAM pair potential 
for several elements. We show values of A, and 5 for 
elements commonly found in BMGs in Table ITTll Pb, Pd, 
Pt, Mg, Fe, Ta, Au, Ti, Mo, W, and Nb do not have 
significant intermediate-range repulsive interactions. 

To study the effects of intermediate-range repulsive in¬ 
teractions on the GFA, we performed MD simulations of 
N = 1372 spherical atoms that interact pairwise via the 
DZ potential. We followed the same cooling protocol as 
used for the simulations of Lennard-Jones systems with 
pressure that decreases exponentially in time as discussed 
in Sec. 12 2.11 We fixed the strength of the intermediate- 
range repulsive interactions at ^ = 0.35e and varied A and 
S to tune the location of the peak {X + 6)/2 and range 
i5 — A of Mbump. We also studied binary mixtures com¬ 
posed of A atoms that interact via the DZ potential with 
^ = 0.35e, A = 1.2, and S = 2.15, and B atoms that in¬ 
teract via the LJ potential with diameter ratio a = 0.95. 
The number fraction of small atoms xb is varied from 0 
to 1 in steps of 0.2. 


2.4. LJ-EAM and EAM potential 

The total potential energy U employed in the 
embedded-atom model for metals includes pairwise and 


TABLE IV: The initial and final temperatures, Ti and Tf, em¬ 
ployed during the cooling protocol in the molecular dynamics 
simulations of the LJ-EAM potential with many-body interac¬ 
tion strength A and electron density inverse decay length /3. 


A (eV) p (A T{K) Tf{K) A (eV) p (A Ti{K) Tf{K) 


0 

4 

2000 

300 

0.66 

2 

2000 

300 

1.32 

2 

2305 

343 

1.98 

2 

3285 

479 

0.66 

4 

2000 

300 

1.32 

4 

2257 

336 

1.98 

4 

3253 

475 

0.66 

6 

2000 

300 

1.32 

6 

2242 

337 

1.98 

6 

3271 

472 


many-body contributions: 

= + ( 8 ) 

i<j i 

where the many-body embedding function Fi depends 
on the electron density associated with each atom i (nor¬ 
malized by el(j\) and pf = ^ HE HI) HE- To 

quantify the effects of the many-body interactions on 
the GFA, we focused on the LJ-EAM potential, where 
uivij) = ULj{rij), Fi{pf) = Ap®(lnp® - rmlcrA)l‘2 and 
p^{rij) = Cexp[—— Vm)], where C and are cal¬ 
ibrated to experimental data on alloys [H, [s^. We set 
the atomic diameter a a = 2.8 A and attraction depth 
e = 0.2 eV for the LJ potential to match the pair po¬ 
tential of typical metals such as Zr. The parameters A 
and j3 control the many-body interaction strength and 
inverse decay length of the electron density, respectively. 

We performed MD simulations of the LJ-EAM for sev¬ 
eral pure metals and of the full EAM for several binary 
alloys using the LAMMPS simulation software [b^. We 
cooled systems in the liquid state to low temperature at 
constant zero pressure at different rates R. The initial 
and final temperatures for several systems (specified by 
A and (3) are summarized in Table IIVI For our studies 
of the full EAM potential, we set N = 4000 and fixed 
the initial and final temperatures at Ti = 2000Ar and 
Tf = 300A:. 


2.5. Critical cooling rate 

To calculate the critical cooling rate Rc for each metal¬ 
lic system, we initialized the liquid state at high tempera¬ 
ture, cooled the system exponentially to low temperature 
at a given rate R at either fixed volume or exponentially 
decaying pressure as in Eq. [31 and measured the global 
bond orientational order parameter Qq Q. For hard- 
sphere interactions, we compressed the systems so that 
the packing fraction approached that at jamming onset 
exponentially, which is thermodynamically equivalent to 
cooling systems exponentially . For all systems stud¬ 
ied, the average global bond orientational order param¬ 
eter Qe versus log R possesses a sigmoidal shape with a 
midpoint defined by Rc- Below, we show results for R^ 
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FIG. 3: The critical cooling rate Rc (in units of 10^® K/s) as a 
function of the repulsive core softness exponent m in monodis- 
perse systems with generalized repulsive Lennard-Jones (tri¬ 
angles) and Lennard-Jones (squares) interactions cooled at 
constant density pa\ = 1 and Lennard-Jones interactions 
with an exponentially decaying pressure given in Eq. [3] (cir¬ 
cles). Variations in the softness exponent lead to different 
crystalline structures that compete with glass formation in¬ 
cluding face-centered cubic (FCC; empty symbols) and body- 
centered cubic (BCC; filled symbols). 


for the pair potentials described in Secs. 12 2.1112 2.3l and 
the full and LJ-EAM potential in Sec. 12 2.41 


3. RESULTS 

3.1. Core Softness 

To investigate the effects of softness of the repulsive 
core on the GFA, we first measured the critical cooling 
rate Rc for monodisperse systems that interact via the 
generalized L J (Eq. [T|) and RLJ (Eq. [2]) pairwise poten¬ 
tials as a function of the softness exponent for m = 1, 
3, 5, 8, 10, and 12. As shown in Eig. [31 when cooling 
at constant number density pa\ = 1, the GEA increases 
weakly [Rc decreases by less than an order of magni¬ 
tude) as the repulsive core becomes softer (to decreases). 
When cooling a LJ system with a pressure that decays 
exponentially in time as in Eq. [31 the dependence of Rc 
on the softness exponent to is even weaker, except for 
systems with extremely soft core repulsions with to = 1. 
In contrast, most atomic species that are found in BMGs 
possess TO > 4 (Table m. 

As shown in Fig.jSl the crystalline structures that com¬ 
pete with glass formation in systems with core-softened 
RLJ interactions at p(j\ = 1 are face-centered cubic 
(FCC) for all exponents to studied. In addition, FCC 
crystals compete with glass formation in LJ systems, but 
as the repulsive core softens, body-centered cubic (BCC) 
crystals become more stable [1^. We find that BCC is 


8 -'- 1 -'- 1 -'- 1 -'- 1 - 1 -r 



0 -^^- 1 -^-'-■-'-^- 1 -- 

0 1 2 3 4 5 6 

Ti/O 


FIG. 4: Radial distribution functions g[rij) (vertically shifted 
for visualization) for monodisperse spherical atoms with di¬ 
ameter o that interact via the generalized m-6 LJ potential 
(Eq.[T]) with m = 1 (top), 3 (middle), and 12 (bottom) cooled 
at rate 7? = 0.1 > Rc- Compared to the m = 12 LJ system, 
the soft to = 1 LJ system shows strong volume contraction 
with the first peak shifted to smaller separations Vij. Systems 
with intermediate softness m — 3 display two isostructural 
states: contracted high density (dotted line) and expanded 
low density (dashed line) glasses. The left and right insets 
show snapshots of the contracted and expanded to = 3 LJ 
systems, respectively. 


the crystal type that competes with glass formation for 
TO = 3 LJ systems cooled at constant density pa\ = 1 
and for to = 3 and 5 LJ systems cooled such that the 
pressure obeys Eq. [31 

Structural characterizations of atomic systems that 
interact via the generalized LJ potential are shown in 
Fig. m for cooling rates R > Rc- As the repulsive core 
of the potential becomes softer (Le. to decreases), the 
attractive well of the potential widens to include second- 
neighbor attractive interactions, which can compensate 
repulsive first-neighbor interactions. Indeed, LJ systems 
with TO = I and 3 exhibit phase separation into dilute and 
compressed regions when cooled at fixed density pa\ = 1 
and volume contraction, where the first neighbor separa¬ 
tions are smaller than the location of the potential min¬ 
imum, when cooled such that the pressure obeys Eq. (3] 
In fact, the to = 3 LJ system displays two isostruc¬ 
tural glassy states, contracted and expanded, with dif¬ 
ferent densities as shown in the inset to Fig. 01 Similar 
isostructural transitions have been found in equilibrium 
systems with narrow-ranged attractive interactions [h^. 
Large density differences between polymorphs in metal¬ 
lic glasses such as those found in CessAUs are often at¬ 
tributed to electronic many-body interactions . How¬ 
ever, here we show that softening the pairwise repulsive 
core (which increases the range of the attractive well) can 
also give rise to polymorphs with different densities. 

We also investigated the effects of core softness on the 
glass-forming ability in binary mixtures that interact via 
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the generalized m-6 LJ potential. We focused on three 
mixtures with diameter ratio a = <JbI<Ja = 0.95: (1) 
conventional LJ systems with m = 12, (2) core softened 
LJ systems with m = 5, and (3) mixtures of LJ systems 
with m = 12 {A species) and m = 5 {B species). While 
FCC is the crystalline structure that competes with glass 
formation for binary LJ systems with m = 12, BCC is the 
competing crystalline structure for binary mixtures with 
TO = 5 for all number fractions Xb cls shown in Fig. [5] For 
both TO = 12 and to = 5 systems, the variation in Rc{xb), 
which is less than an order of magnitude, is controlled by 
the diameter ratio a = 0.95. In binary mixtures of LJ 
systems with to = 12 and to = 5 interactions, FCC re¬ 
mains the crystalline structure that competes with glass 
formation, except when xb ~ 1. However, because of 
the incompatibility between FCC and BCC crystalline 
structures, the GFA for the m = 12 and to = 5 LJ mix¬ 
tures is significantly enhanced compared to glasses with 
TO = 12 or TO = 5 interactions alone. For example, Ni-Ta 
is a good glass former despite the fact that it possesses 
a diameter ratio near unity (a ~ 0.9) (^ . Incompatibil¬ 
ity between competing BCC and FCC crystal structures 
is a possible cause of the enhanced GFA. As shown in 
Table m Ni has a relatively large pairwise repulsive ex¬ 
ponent (6 < TO < 10) with equilibrium FCC structure, 
while Ta has a relatively small exponent (3 < to < 5) 
with equilibrium BCC structure Since the softness 
exponents of the pairwise interactions vary significantly 
from one element to another (Table [J), softness-induced 
competing crystal incompatibility can enhance the GFA 
of binary and multi-component BMG-forming alloys. 


3.2. Non-additivity 

We performed event-driven molecular dynamics sim¬ 
ulations of binary non-additive hard spheres (Sec. 12 2.21) 
to investigate the effects of non-additivity of the pairwise 
repulsive interactions on the GFA of alloys. We measured 
the critical cooling rate Rc of non-additive binary hard 
spheres with diameter ratios a = gbIcfa = 1-0, 0.97, 
0.95, 0.93, 0.9, and 0.5 and number fractions of the small 
spheres xb = 0.5 and 2/3 over a range of non-additivity 
parameters S. Since S > 0 is rare among binary alloys 
f Table Hill, we expect that hard-sphere systems with pos¬ 
itive non-additivity are poor glass-formers. For example, 
we find that systems with a = 1 and S = 0.05 display 
strong demixing between A and B particles and are not 
good glass formers. 

Our previous studies of additive binary hard spheres 
(S = 0) have shown that well-mixed FCC solid solutions 
are the crystal structures that compete with glass forma¬ 
tion when a > 0.8, while the systems tend to demix when 
a < 0.8 [l^. For E < 0 and a = 1.0, 0.97, 0.95, 0.93, and 
0.9, the GFA improves as S becomes more negative, and 
the competing crystal structure remains the FCC solid 
solution (Fig. [5]). The change in Rc with decreasing S 
also increases as a decreases with roughly an order of 



FIG. 5: The critical cooling rate Rc. for binary mixtures of 
spherical atoms that interact via the generalized m-6 LJ po¬ 
tential (Eq.[T]) at diameter ratio obI^xa = 0.95 is plotted as a 
function of the number fraction of small atoms xb- We show 
to = 5 (squares) and 12 (triangles) LJ systems, mixtures (di¬ 
amonds) of TO = 12 (species A) and 5 (species B), as well as 
mixtures (circles) of spheres with DZ (A species) and m = 12 
LJ [B species) interactions. Open and filled symbols indi¬ 
cate that the crystalline structure that competes with glass 
formation is FCC and BCC, respectively. 

magnitude difference in Rc between systems with E = 0 
and E = —0.05 at a = 0.9. Enhancement of the GFA 
arising from non-additivity of the repulsive cores (E < 0) 
has also been observed in LJ systems [^ . 

For binary systems with large atomic size differences 
(j.e. a <C 0.8), the variation of Rc with E is opposite 
to that obtained for binary systems with small atomic 
size differences. As shown in Fig. [SI we find that Rc 
grows with increasing E at a = 0.5. For a = 0.5 and 
E < 0, compound crystals are the ordered structures 
that compete with glass formation since negative non¬ 
additivity promotes mixing. As an example, although 
the AB 2 compound is the densest crystal for binary hard 
spheres with a = 0.5 and E = 0, it is not kinetically 
accessible during compression due to the strong drive for 
demixing [H, [^, . However, when E becomes nega¬ 

tive {e.g. E = —0.05), we find that the AB 2 compound 
forms easily for the compression rates that we studied, 
as shown in the inset to Fig. jS] Thus, the formation of 
intermetallic compounds in alloys can be enhanced by 
pairwise negative non-additivity among different atomic 
species. 

3.3. Intermediate-range Repulsive Interactions 

We also investigated crystallization and glass forma¬ 
tion as a function of the form of intermediate-range repul¬ 
sive pairwise interactions fSec. 12 2.^ . We first performed 
molecular dynamics simulations of monodisperse spheres 
interacting via the DZ potential (Eq. (5]) at fixed strength 
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FIG. 6: Critical cooling rate Rc (in units of 10^^ K/s) plot¬ 
ted as a function of the nonadditivity parameter E for bi¬ 
nary hard-sphere systems for several diameter ratio and small- 
sphere number fraction combinations: a = 1.0 and xb = 0.5 
(crosses), 0.97 and 0.5 (squares), 0.95 and 0.5 (pentagons), 
0.93 and 0.5 (downward triangles), 0.90 and 0.5 (diamonds), 
and 0.5 and 2/3 (upward triangles). The inset shows snap¬ 
shots of (top) demixed and (b) compound crystals that form 
for i? < i?c at E = 0.0 and —0.05, respectively, with diameter 
ratio a = 0.5. 


^ = 0.35e and varying peak location (A + 5)/2 and width 
(5 - A. In Fig. [3 we plot the critical cooling rate Rc as a 
contour plot versus {\+5)/2 and 5—X over ranges that are 
relevant to BMGs (Table Hill) . We find several regions of 
good glass-forming ability (small Rc) and different crys¬ 
tal structures that compete with glass formation. For a 
large region of parameter space, FCC is the competing 
crystal structure. BCG is the competing crystal structure 
when the location of the peak in Ubump approaches third- 
neighbor separations at Ri Virm- We also find an “8- 
4” crystal structure that competes with glass formation, 
with atom positions located on embedded octagons and 
squares when they are projected into two dimensions. 
(See the inset of Fig. (T]). In three dimensions, one can 
see that the atoms forming the octagons and squares are 
located in alternating stacked layers. (See Fig. [5] for a 
comparison of the radial distribution functions for FCC, 
BCG, and 8-4 crystals.) When the intermediate-range 
repulsion becomes too strong (ie. large S), microphase 
separation becomes energetically favorable compared to 
macroscale phase separation [^, . 

We also studied the critical cooling rate Rc for bi¬ 
nary mixtures {e.g. Zr-Cu alloys), in which one compo¬ 
nent possesses intermediate-range repulsive interactions 
and the other component does not. We focused on bi¬ 
nary systems with atoms that interact via the DZ [A 
species) and LJ potential [B species) with diameter ratio 
(^bI^a = 0.95. For the DZ potential, we set the param¬ 
eters C/£ ~ 0.4, A ~ 1.2, and 8 ~ 2.2 to mimic those 
of Zr atoms (Table Hill) . As shown in Fig. O Rc for this 
binary mixture is suppressed by more than two orders 
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FIG. 7: Gontour plot of the critical cooling rate Rc (in reduced 
unit) for monodisperse spheres that interact via the DZ poten¬ 
tial (Eq. [6| as a function of the location of the peak (A -|- (5)/2 
and width J — A of Ubump- The bounds for the parameters are 
determined by A > rmjoA = 1.12 and J < rdcA = 2.5. Rc 
contours are interpolated from simulation data points. The 
symbols indicate where FGC (triangles), BGC (squares), 8-4 
(stars) crystalline structures, and microphase separation (cir¬ 
cles) is observed. Crosses indicate systems for which the com¬ 
peting crystal structure is unknown and Rc is estimated from 
the slowest cooling rate employed. The inset shows a snapshot 
of a 8-4 crystal that includes top (dark) and bottom (light) 
layers of atoms with square symmetry (red squares). 


of magnitude compared to the pure system with LJ or 
DZ interactions alone because the two species possess in¬ 
compatible equilibrium crystal structures {i.e. FCC and 
BCG). This mechanism of incompatible equilibrium crys¬ 
tal structures may explain the exceptionally good glass¬ 
forming ability of the Zr-Cu system, even though it is a 
binary, rather than, multi-component alloy. 


3.4. LJ EAM for Monoatomic Systems 

To determine the relative contributions of the pairwise 
and many-body interactions to the GFA of alloys, we per¬ 
formed molecular dynamics simulations of the LJ-EAM 
potential (Sec. 12 2.41) as a function of the many-body in¬ 
teraction strength A and electron density inverse decay 
length P for monoatomic systems. In Fig. 13 we show the 
critical cooling rate Rc for monodisperse LJ-EAM sys¬ 
tems as a function of A for /3 = 2, 4, and 6 A . We find 
that Rc ~ 10^^ K/s. Rc changes by less than one order 
of magnitude as A and /3 are varied over the range that 
is relevant for elements found in BMGs even though the 
total potential energy per atom U/N varies linearly with 
A. We also find that FCC crystals are the ordered struc¬ 
tures that compete with glass formation in monoatomic 
LJ-EAM systems over the full parameter range for A and 
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FIG. 8: Radial distribution function g(rij) (solid lines, left 
axis) and pair potential u(rij) (dashed lines, right axis) for 
monodisperse spheres that interact via (a) the Lennard-Jones 
potential in a FCC crystal structure and via the DZ potential 
with (b) A = \/2 and S — 2.1 in a FCC crystal structure, (c) 
A = 1.2 and <5 = 2.3 in a BCC crystal structure, and (d) A = 
1.5 and <5 = 2.2 in a 8-4 crystal strncture. The vertical dotted 
and dashed lines indicate the BCC lattice spacings (relative 
to Tm) l:2/'\/3:2 and the FCC lattice spacings l:%/2:%/3:2 up 
to third and fourth nearest neighbors, respectively. 


/3. Thus, we argue that many-body interactions have a 
weak influence on the GFA compared to the pairwise in¬ 
teractions for inonoatomic systems. 


3.5. Full EAM for Binary Alloys 

We also measured the critical cooling rate Rc for sev¬ 
eral binary alloys as a function of the number fraction 
xb of the small atomic species using the full EAM po¬ 
tential. We focused on Zr-Cu, Mg-Al, and Cu-Nl alloys 
with atomic diameter ratios that range from a = 0.79 to 
0.98. In Fig. 1101 we compare Rc versus Xb from simu¬ 
lations of the full EAM potential for these alloys to Rc 
obtained from simulations of additive hard spheres with 
comparable values of a [T^. 

As expected, Rc for binary alloys with a ~ 1 {i.e. 
Cu-Nl) Is nearly Independent of xb- In addition, when 



FIG. 9: The critical cooling rate Rc from simulations of the 
LJ-EAM plotted as a function of the many-body interaction 
strength A (in eV) for several values of the electron density 
inverse decay length /3 = 2 (squares), 4 (circles), and 6 A 
(triangles). Rc from simulations of the full EAM for Zr (i.e. 
A « 1.32 eV and /3 ~ 4 A ) is indicated by the horizontal 
dashed line. Error bars give the standard deviation from 10 
independent simulations with random initial conditions. The 
inset shows the total potential energy U/N per atom (for 
cooling rates R > Rc) versus A for the same values of /3 in 
the main panel. U/N for the full EAM of Zr is given by the 
horizontal dashed line. 


the hard-sphere simulations with a = 1 are calibrated 
to Ni, Rc from simulations of the hard-sphere and EAM 
potentials agree semi-quantitatively. From our previous 
simulations of hard spheres [l^, we know that Rc{xb) 
develops a deep minimum that shifts to larger xb ss a 
decreases from unity. For example, when a = 0.9, Rc 
for hard-sphere systems at xb ~ 0.6 is two orders of 
magnitude less than the value when a = 1. Although we 
are not able to simulate sufficiently slow rates, it appears 
that Rc at the minimum in xb for Mg-Al with a = 0.94 
will decrease by at least two orders of magnitude and the 
minimum in Rc{xb) will occur at xb > 0.5. We also find 
similar results for Rc for hard spheres with a = 0.79 and 
for EAM of Zr-Cu with a deep minimum in the range 
0.2 <XB < 0 . 8 . 

We also determined the crystal structures that com¬ 
pete with glass formation in the full EAM simulations 
of binary alloys. We find that FCC (or HCP) is most 
often the competing crystal structure, as in simulations 
of additive binary hard spheres, but we also find excep¬ 
tions. In particular, we show that on the Zr-rich side of 
Zr-Cu, BCC crystal structures compete with glass forma¬ 
tion. The BCC equilibrium structure for the Zr-Cu alloys 
can likely be attributed to the pairwise part of the EAM 
potential. For example, the pair potential for Zr pos¬ 
sesses intermediate-range repulsive interactions with the 
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FIG. 10: The critical cooling rate Rc (in K/s) for several 
binary alloys, including Zr-Cu with atomic diameter ratio 
a = 0.79 (circles), Mg-Al with a = 0.94 (triangles), and Cu-Ni 
with a = 0.98 (diamonds), is plotted as a function of number 
fraction of small atoms xb using molecular dynamics simu¬ 
lations of the full EAM. Error bars on Rc are obtained from 
the standard deviation from 5 independent simulations. The 
EAM source files are given in Refs. |^ . [^ . . As a com¬ 

parison, Rc for additive binary hard spheres with a = 1.0, 
0.9, and 0.79 are shown as dashed lines. We also indicate 
when FCC or HOP (open symbols) and BCC (filled symbols) 
crystal structures compete with glass formation. 


location of peak (A-b :5)/2 = 1.70 and width 5 — \ = 1.08 
(Table IIIII) in a region of parameter space that has been 
shown to display BCC crystal structure (Fig. [T]). 


4. CONCLUSION 

The hard-sphere model has provided a predictive de¬ 
scription of cry stallization and glass formation in sim¬ 
ple liquids [33|. In addition, we have shown in recent 
studies that the additive hard-sphere model can explain 
more than 13 orders of magnitude variation in the critical 
cooling rate Rc, which nearly spans the full range of GFA 
from that for pure metals to that for the best BMGs [l^ . 
We also showed that the best binary and ternary BMGs 
occur in the region of parameter space (i.e. diameter ra¬ 
tio and number fraction) with the smallest values of Rc 
for hard spheres. 

However, in metallic systems, there are a number of 
additional features of the interatomic potential beyond 
hard-core repulsions, including softness, non-additivity, 
and range of the pairwise interactions. For example, 
metallic atoms typically appear softer (with smaller val¬ 
ues of the exponent of the repulsive core) than the com¬ 
monly used LJ pair potential and possess several per 


cent negative non-additivity due to shortening of metal¬ 
lic bonds [ 13 . In addition, Friedel oscillations in metals 
give rise to intermediate-range repulsion at separations 
beyond the short-range attractive well [s^. The inter¬ 
atomic potential for metals also includes many-body in¬ 
teractions from the electronic degrees of freedom. In this 
manuscript, we investigated how these additional fea¬ 
tures affect the GFA of pure and binary metallic systems. 

We performed molecular dynamics simulations of sev¬ 
eral model systems to study the effects on the GFA 
for each of the key features of the interatomic poten¬ 
tial separately. For example, we performed simulations 
of monodisperse and binary spheres that interact via the 
generalized LJ and DZ pair potentials to quantify the ef¬ 
fect of the softness of the repulsive core and form of the 
intermediate-range repulsive interactions on the GFA. 
We also performed MD simulations of non-additive bi¬ 
nary hard spheres to quantify the effects of non-additivity 
on the GFA. We found that softness, non-additivity, and 
form of the intermediate-range repulsions cause devia¬ 
tions in Rc that are only 1 ~ 2 orders of magnitude from 
the additive hard-sphere predictions. 

While FGG is the most stable crystal structure for LJ 
and hard-sphere systems, softening of the repulsive core 
gives rise to novel contracted disordered structures, as 
well as the formation of BGG crystals. We also showed 
that negative non-additivity of the repulsive core in bi¬ 
nary alloys improves the GFA when the competing crys¬ 
tal structures are solid solutions. However, when the 
atomic size ratio is in the demixing regime (a < 0.8), neg¬ 
ative non-additivity can favor the formation of compound 
crystals and decrease the GFA. The crystal structure that 
competes with glass formation, and thus the GFA, also 
depends sensitively on the form of the intermediate-range 
repulsive interactions. We find that when the competing 
crystal structures of each component in an alloy are in¬ 
compatible {e.g. FCC and BCC), the GFA can be en¬ 
hanced compared to hard-sphere predictions. 

We also investigated the relative contributions of the 
pairwise and many-body interactions to the GFA by per¬ 
forming molecular dynamics simulations of the LJ-EAM 
potential. We found that including the many-body inter¬ 
actions only changes Rc by less than one order of mag¬ 
nitude compared to that when the many-body interac¬ 
tions are not included. We also calculated Rc for several 
binary alloys using the full EAM potential and found 
qualitatively the same results as for binary hard spheres. 
Thus, we argue that hard-sphere interactions provide a 
qualitatively accurate model for predicting the GFA of 
alloys. Other features of the interatomic potential (be¬ 
yond additive hard-core repulsion) give rise to only 1-2 
orders of magnitude variation of Rc, which is small com¬ 
pared to the more than 13 orders of magnitude variation 
predicted by hard-sphere systems. Despite this, includ¬ 
ing additional features to the interatomic potential be¬ 
yond hard-sphere interactions is important for the design 
of new BMGs since precise quantification of the critical 
casting thickness can determine whether a new BMG is 
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commercially viable. 
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